

clear
capture log close
set more off

global path ""
cd ""
global dta1 ""


clear
use income_b1990.dta // after transformed income

forvalues i=1990/2020 {
use income_b1990.dta

* create real income groups

	gen hhfaminc_1990 = 1 if hhincome_1990<3000
	replace hhfaminc_1990 = 2 if hhincome_1990>=3000 & hhincome_1990<4000
	replace hhfaminc_1990 = 3 if hhincome_1990>=4000 & hhincome_1990<5000
	replace hhfaminc_1990 = 4 if hhincome_1990>=5000 & hhincome_1990<6000
	replace hhfaminc_1990 = 5 if hhincome_1990>=6000 & hhincome_1990<7500
	replace hhfaminc_1990 = 7 if hhincome_1990>=7500 & hhincome_1990<9000
	replace hhfaminc_1990 = 8 if hhincome_1990>=9000 & hhincome_1990<10000
	replace hhfaminc_1990 = 9 if hhincome_1990>=10000 & hhincome_1990<11000
	replace hhfaminc_1990 = 10 if hhincome_1990>=11000 & hhincome_1990<12500
	replace hhfaminc_1990 = 12 if hhincome_1990>=12500 & hhincome_1990<14000
	replace hhfaminc_1990 = 13 if hhincome_1990>=14000 & hhincome_1990<15000
	replace hhfaminc_1990 = 14 if hhincome_1990>=15000 & hhincome_1990<17500
	replace hhfaminc_1990 = 15 if hhincome_1990>=17500 & hhincome_1990<20000
	replace hhfaminc_1990 = 16 if hhincome_1990>=20000 & hhincome_1990<22500
	replace hhfaminc_1990 = 17 if hhincome_1990>=22500 & hhincome_1990<25000
	replace hhfaminc_1990 = 18 if hhincome_1990>=25000 & hhincome_1990<27500
	replace hhfaminc_1990 = 19 if hhincome_1990>=27500 & hhincome_1990<30000
	replace hhfaminc_1990 = 20 if hhincome_1990>=30000 & hhincome_1990<32500
	replace hhfaminc_1990 = 21  if hhincome_1990>=32500 & hhincome_1990<35000
	replace hhfaminc_1990 = 22 if hhincome_1990>=35000 & hhincome_1990<40000
	replace hhfaminc_1990 = 23 if hhincome_1990>=40000 & hhincome_1990<50000
	replace hhfaminc_1990 = 24 if hhincome_1990>=50000 & hhincome_1990<75000
	replace hhfaminc_1990 = 25 if hhincome_1990>=75000 & hhincome_1990!=.

keep if year==`i'

collapse (sum) asecfwt, by(hhfaminc_1990)
rename asecfwt hh_count

egen sumhh=total(hh_count)
gen weight=hh_count/sumhh
save `i'_cps_income,replace
}

*** predictefd kwh

clear
forval i=1991/2020  {
clear
use recs_1990
capture drop if kwh==999999

quietly bysort moneypy: egen sumhh=total(nweight)
quietly bysort moneypy: gen weight=nweight/sumhh
quietly bysort moneypy: egen meankwh=total(weight*kwh)

collapse meankwh,by(moneypy)

quietly drop if moneypy==. 
rename moneypy hhfaminc_1990
sort hhfaminc_1990 
quietly destring hhfaminc_1990,replace

quietly merge 1:1 hhfaminc_1990 using  `i'_cps_income

quietly drop if meankwh==.
quietly drop if hh_count==.

*sum kwh // this sum is not estimated/ estimated is by weights, this is just summary
quietly sum meankwh[aweight=weight]
quietly return list
dis `i',r(mean)

}


***standard errors

*****====standard errors====*****

clear
use recs_1990

sort moneypy division

egen hhgroup=group(moneypy)

reg kwh i.hhgroup [weight=nweight], nocons 
estimates store results1
esttab results1, b(3) se(3) scalars(r2) star( * 0.10 ** 0.05 *** 0.01) 


matrix B=e(V)

local dim (`= rowsof(A)',`=colsof(A)') 

di "`dim'" 

forval i=1990/2020 {
quietly use `i'_cps_income,clear
quietly mkmat weight,matrix(A)
quietly matrix SD=A'*B*A
quietly matrix list SD
quietly scalar z=SD[1,1]
di sqrt(z)
}


**** predicted growth


**** generate cps_division count

clear
use income_b1990.dta
bysort year: egen sumhh=total(asecfwt)
bysort year: gen weight=asecfwt/sumhh
bysort year: egen meanincome=total(weight*hhincome_1990)

collapse meanincome,by(year)

gen grate=meanincome/meanincome[_n-1]
replace grate=1 if missing(grate)



gen lngrate = ln(grate)
gen lntotal = sum(lngrate)
gen double grateproduct = exp(lntotal)
drop grate
rename grateproduct grate
drop lngrate
drop lntotal
drop meanincome

save grate_1990.dta,replace  //FIGURE 3


clear
use income_b1990.dta
merge m:1 year using grate_1990 
drop _merge
sort year

gen nn_income=hhincome_1990*grate
order hhincome nn_income grate

save income_b1990_grate.dta,replace
** 1990 -2020
clear

forvalues num=1990/2020 {
use income_b1990_grate.dta
	*** new income groups
	
		gen ngroup = 1 if nn_income<3000
	replace ngroup = 2 if nn_income>=3000 & nn_income<4000
	replace ngroup = 3 if nn_income>=4000 & nn_income<5000
	replace ngroup = 4 if nn_income>=5000 & nn_income<6000
	replace ngroup = 5 if nn_income>=6000 & nn_income<7500
	replace ngroup = 7 if nn_income>=7500 & nn_income<9000
	replace ngroup = 8 if nn_income>=9000 & nn_income<10000
	replace ngroup = 9 if nn_income>=10000 & nn_income<11000
	replace ngroup = 10 if nn_income>=11000 & nn_income<12500
	replace ngroup = 12 if nn_income>=12500 & nn_income<14000
	replace ngroup = 13 if nn_income>=14000 & nn_income<15000
	replace ngroup = 14 if nn_income>=15000 & nn_income<17500
	replace ngroup = 15 if nn_income>=17500 & nn_income<20000
	replace ngroup = 16 if nn_income>=20000 & nn_income<22500
	replace ngroup = 17 if nn_income>=22500 & nn_income<25000
	replace ngroup = 18 if nn_income>=25000 & nn_income<27500
	replace ngroup = 19 if nn_income>=27500 & nn_income<30000
	replace ngroup = 20 if nn_income>=30000 & nn_income<32500
	replace ngroup = 21  if nn_income>=32500 & nn_income<35000
	replace ngroup = 22 if nn_income>=35000 & nn_income<40000
	replace ngroup = 23 if nn_income>=40000 & nn_income<50000
	replace ngroup = 24 if nn_income>=50000 & nn_income<75000
	replace ngroup = 25 if nn_income>=75000 & nn_income!=.



keep if year==`num'

collapse (sum) asecwth, by(ngroup)
rename  asecwth hh_count

egen sumhh=total(hh_count)
gen weight=hh_count/sumhh
save grate`num'_cps_90,replace
}


*** predictefd kwh

clear
forval i=1990/2020  {
clear
use recs_1990 //
** get weighted kwh

drop if moneypy==.
bysort moneypy: egen sumhh=total(nweight)
bysort moneypy: gen weight=nweight/sumhh
bysort moneypy: egen meankwh=total(weight*kwh)

collapse meankwh,by(moneypy)

rename moneypy ngroup
quietly sort ngroup 
quietly destring ngroup,replace

quietly merge 1:1 ngroup using  grate`i'_cps_90

quietly drop if meankwh==.
quietly drop if hh_count==.

quietly sum meankwh[aweight=weight]
quietly return list
dis `i',r(mean)

}


*** standard errors of predicted growth
clear
use recs_1990

sort moneypy division

egen hhgroup=group(moneypy )

reg kwh i.hhgroup [weight=nweight], nocons 
estimates store results1
esttab results1, b(3) se(3) scalars(r2) star( * 0.10 ** 0.05 *** 0.01) 


matrix B=e(V)

local dim (`= rowsof(A)',`=colsof(A)') 

di "`dim'" 

local dim (`= rowsof(B)',`=colsof(B)') 

di "`dim'" 

forval i=1990/2020 {
quietly use grate`i'_cps_90,clear
quietly mkmat weight,matrix(A)
quietly matrix SD=A'*B*A
quietly matrix list SD
quietly scalar z=SD[1,1]
di sqrt(z)
}
